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Abstract: The field of proton accelerators research is experiencing a surge of interest as proton beams find 
wide applications in areas like cancer therapy, industrial applications, and scientific researches. The recent 
advancements in high gradient radio frequency standing wave and travelling wave ion accelerating structures, 
reaching frequencies of several GHz in the past few years, have paved the way for compact linear ion 
accelerators with high gradients. Compared with the abundant design code for standing wave ion linacs, the 
design code for travelling wave structure is very limited. We developed a new type of particle dynamics 
tracking (PDT) code, which can be used in beam dynamics design for both standing wave and travelling wave 
structures based on the electromagnetic field distribution of the unit cell. The benchmark with the widely used 
commercial linac beam dynamics simulation code TraceWin shows its accuracies both in energy gain and 
envelopes are good. The code has been applied in the beam dynamics design of a low current backward 
travelling wave proton linac, which operates at RF frequency of 2.856 GHz and working at 52/6 mode. The 
linac spans a total length of approximately 3 m, comprising 10 tanks and 10 PMQs (Permanent Magnet 
Quadrupoles). It operates with a 30 MV/m accelerating gradient, characterized by full transmission and minimal 
emittance growth. 
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1 Introduction 


Way back in the beginning of the last century, the proton accelerator has been put forward and operated successfully, 
and various types of accelerators have been clearly thriving ever since. The accelerated particle beams have a wide range 
of applications and prospects in many areas such as nuclear physics [1], materials science, space exploration and 
especially in radiation treatment [2-3]. In recent years, proton and heavy ion radiotherapy have seen rapid development 
worldwide. At the time of writing, January 2023, there are around 120 particle therapy facilities in operation and most of 
them are based on cyclotrons and synchrotrons [4]. 

The frequency of proton linear accelerators (linacs) are typically 200-400 MHz, with a low acceleration gradient 
and large volume, coupled with large-scale investments. This is one of the limiting factors that limits proton therapy from 
becoming the standard of broader cancer care at present. High accelerating gradient ion linacs have become feasible 
solutions with S-band frequencies, which also provides a new idea for the application of proton linac. A research group 
designed the LIBO (LInac BOoster), a side-coupled Linac that works at 3 GHz and whose effective accelerating gradient 
can reach 28.5 MV/m, can be used downstream of a linac or a cyclotron, to boost the low energy of proton beam from 
about 60 MeV to 200 MeV or more[5-6]. ACCIL is a compact ultra-high gradient coupled cavity linac (CCL) acceleration 
structure that operates at 2856 MHz and can accelerate up to 50 MV/m. Its RF parameters degrade dramatically for 
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structures with £<0.7 [7-8]. The TULIP project is dedicated to the development of a compact proton linac. A new 
backward travelling waveguide (BTW) RF structure with a high gradient is designed, which can meet the accelerating 
gradient of 50 MV/m. The low shunt impedance caused by low-f CCL structure is solved by adding external disc coupling 
holes for magnetic coupling. In order to achieve low peak electric field and high shunt impedance with lower 2, the article 
[9] proposed to synchronize beam with the -1 harmonic wave, which will prolong the acceleration period and increase 
the shunt impedance without significantly increasing the peak field. 

While TraceWin is capable of simulating beam dynamics through fixed standing wave structures, it may not provide 
extensive design support for travelling wave RF structures [10]. In the literature [11-12], a structure optimization by 
polynomial fitting (SOPF) for standing wave structures is compiled, which is a multi-objective optimization algorithm. 
The frequency sensitive values were obtained by the scanning points, and a complete database was established by 
polynomial fitting according to the analysis results, and the relevant parameters were finally determined. There are also 
codes that can track particles through travelling wave structures, such as RF-Track [13-14], which can be used to 
maximize the transmission effect as well as validate the lattice design. 

Compared with the abundant available design tools for standing wave proton linac, the design tool for tracking 
protons through travelling wave structures is relatively hard to find. Since the electromagnetic structure design and 
particle dynamics design are closely related, it is important to find a starting point. It is urgently needed to compile the 
particle dynamics tracking (PDT) code of our own, which can meet the needs, including precise tracking of single-particle 
and beams of protons, as well as automatically optimization tuning of cavity structural parameters. Moreover, we can 
add functions to the system according to our specific requirements flexibly, rather than using a stereotyped black-box 
tool. It can perform dynamics simulation for not only a single acceleration structure, but also for different situations under 
different types of acceleration structures, modes and frequencies, etc. At the same time, a proton linac based on high 
gradient backward travelling wave is designed by using this new tracking code. This linac requires minimum RF design 
input (cell optimization results), and its application on the 100 MeV BTW proton linac beam dynamics simulation is 
presented. 

Except the introduction, the paper is divided into three parts: Firstly, the dynamic equations and algorithms outline 
used in PDT code are introduced in detail, including the function expression of travelling wave field in periodically 
loaded waveguide and its establishment with the help of CST, longitudinal and transverse particle tracking and beam 
matching, as well as the automatic tuning of RF waveguide structure. In the second part, we make a benchmark of PDT 
code with TraceWin, and the simulation results are in good agreement, which verifies the feasibility of PDT. The third 
part is to simply design a multi-tanks linac with an acceleration gradient of 30 MV/m to boost the energy of 30 MeV up 
to 102.9 MeV in the case of protons, with a RF frequency of 2.856 GHz and phase advance of 52/6. A preliminary 


dynamics study has been carried out using PDT. 


2 Particle dynamics tracking and optimization by PDT code 


To study the dynamic process of protons in the travelling wave accelerating structures, it is necessary to establish a 
physical model of the electromagnetic field first. The typical modeling method is to build a three-dimensional cavity 
model, and the widespread availability of commonly used simulation software such as SUPERFISH [15] and CST studio 
suite [16] makes it feasible to extend the use of these programs to travelling wave structures. 

The whole solving process is mainly divided into two modules. The simulation model integrative framework is 
proposed as in Fig.1. The first step is to make full use of automated MATLAB [17] and CST co-simulation to track the 
dynamics of a single particle in the longitudinal direction and finish the main accelerating cavity design (RF structure 
tunning). The linac typically comprises accelerating tanks, and each tank is composed of multiple cells of equal length. 
The particle velocity, 6., within each tank remains almost constant and changes very little, and the corresponding cell 
length change is less than 0.1 mm. However, the J; value varies significantly between tanks, and this cannot be 


disregarded. As a result, the length of the cells in the adjacent tank needs to be adjusted based on the £; value. We utilize 


CST to model and tune the single cell with periodic condition calculating the electromagnetic field on the central axis. 
Structure parameters such as iris aperture radius, cavity radius and the radial distance of coupling holes, etc. are calculated 
by CST optimizer combined with the required RF characteristics and cell length. After calculating the electric field 
distribution of the unit cell in first tank with specific geometric parameters, dynamics tracking of the particles through 
the tank is carried out in MATLAB. Compute the structure parameters of the unit cell in the next tank, such as the cell 
length, tube length even the optimal cell numbers through the particle output state parameters in end of last tank, and take 
them as new input parameters of CST simulation. This cycle continues until the particles are accelerated to the specified 
energy. 

On account of the multi-parameter optimization of the unit cell in each tank must be carried out strictly step by step 
manually adjusting, making the whole process very slow. The tuning work and stored procedure involved can be very 
complex and time-consuming. In order to enable the dynamic invocation, automating RF structure geometric parameters, 
the code has added this module to create appropriate wires to enable communication between PDT and CST, for speeding 
up the design process. 

Into the second module, beam dynamic design is carried out. Including transverse focusing of FODO-like lattice 
approximate matching problem, and finally obtaining the optimal Twiss parameters and lens gradients to achieve 


maximum transmission. As well as longitudinal and transverse multi-particle tracking. 
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Fig.1 (color online) The Simplified scheme of PDT code architecture 
2.1 Travelling wave field in a periodically loaded waveguide 


By CST, the longitudinal field of travelling wave along travelling wave accelerating structure at a certain time can 
be calculated, and then the longitudinal and transverse fields at any other time are represented by transformation methods 
presented in the following sections, which avoids the heavy calculation process and greatly reduces the amount of 
calculation. 

The total field of the lowest order passband of TMO1 mode in periodically loaded waveguide can be expressed as 
the sum of all harmonic fields [18-19], and simplified as: 

E,(r, z, t) = A,(r, z)ei(@t-2@2)) (1) 

The total field is simplified into the exponential form of a complex variable function, that is, the function of the 
amplitude A, and phase 人 办 of the electric field. The field satisfies Maxwell's equations and boundary conditions at any 
time. By obtaining the field along the cavity axis at a specific time to, we can determine the longitudinal field of the 
traveling wave at any other time. For a multi-cell structure, according to Floquet’s theorem, the fields of the cells differ 
only by a constant complex factor, which is denoted as the periodic phase shift ø. Then, we can obtain the field along the 
multi-cell periodic loaded waveguide at any time based on the amplitude A, and phase ¢-(to) = wto-$: function obtained 
from this single cell field calculation. And the field of cell n with equal length / appears as follows: 

E,(z+nl,t + to) = A,e/ (not ot-Pzlto)) (2) 
where E-, A, and pz are all functions that change with z. In the case of a proton machine, the constant gradient mode is 
commonly employed, where the cell length is maintained constant within a short section of the waveguide for a 


compromise between simple fabrication and a little cost of a phase slippage and efficiency loss [20]. There is coupling 


effect between the transverse field and the longitudinal field in the passive case according to Maxwell's equations and 
axial symmetry [21], the transverse electric and magnetic fields under the paraxial approximation can be represented by 
the longitudinal electric field: 


= _OE,r = ae [ee ei(not+wt—gz(to)) Ee ellnorot-extto)| 
i dz 2 2| dz dz < (3) 
r OE, rw 


Bo = =f, te ere) 

As a result, according to the longitudinal electric field at a certain time within a unit cell, the electric field value and 
angular value of magnetic field at any position near the axis at any time in the acceleration gap can be obtained. Based 
on the established travelling wave field, we can further analyze and discuss the transverse and longitudinal dynamic 


equations of particles, and solve the bunched beam problem in periodic focusing channel. 
2.2 Longitudinal dynamic equations and the algorithm 


The velocity vz of the proton passing through the cavity will change with the longitudinal distance z. Let the relative 
phase of the particle entering the cavity be ĝo and the phase of the time component can be expressed as the position 
component az. And the voltage Yoo) felt by the particle through the whole cavity is the integral of the longitudinal electric 
field along the path. Obviously, we can find a 6; that V(0;) reaches a maximum, and particles gain greatest acceleration 
with this incident phase also called synchronize 0° phase. However, in order to accelerate the particles stably and 
considering the particles capture efficiency, it’s necessary to choose a synchronous phase for achieving a compromise 
between less particle loss and higher acceleration efficiency. 

According to the conservation of energy, when the acceleration structure is in TMO1 mode, only the longitudinal 


electric field contributes to the energy gain, and the longitudinal motion equations can be expressed as: 


dy qa, 

a ae? cos (a, — pz + 9) 本 
dwt w 1 

dz c h 


where the independent variable is z, ois the initial input phase. In order to achieve the maximum acceleration efficiency, 
it is crucial for the speed of the travelling wave in the accelerating cavity to match the speed of the particles. This 
relationship is described in detail in section 4. When the number of particles is increased and the initial phase difference 
and energy difference are within a certain range relative to the reference particle, each particle can be tracked and 
positioned to obtain longitudinal phase space diagram of the particles. 

The amplitude and phase of the electric field of unit cell were derived from CST simulation and then imported into 
MATLAB for further analysis. The leapfrog algorithm was employed to solve the longitudinal dynamic differential 
equations of the particles, with the appropriate step threshold selected, to calculate and record the longitudinal velocity 
state quantity f: and the phase state 办 of particles simultaneously. Leapfrog algorithm is fast and accurate enough. And 
cyclic optimization algorithm is used to calculate the synchronize 0° phase 6j, that is, the incident phase when the particle 
gets the maximum energy gain. 

After tracking and recording the state quantities of the longitudinal motion of a single particle, the tracking 
calculation for multi-particles is performed based on the incident condition determined by the reference particle, which 


corresponds to the one with the optimal incident phase 6; that we previously sought. 


2.3 Transverse dynamic equations and algorithm 


The transmission efficiency is closely related to the transverse emittance of the beam. In the context of a linac 


accelerating particles longitudinally, it is important to consider the impact of the cavity on the transverse dynamics of the 


particles. The electromagnetic field can induce defocusing effects on the particles in the transverse direction, potentially 
leading to an increase in the particle's radial distance from the ideal beam path. Excessive radial distances can result in 
particle loss and decreased transmission efficiency. To study the transverse dynamics of a single particle, the transverse 
dynamic equations are analyzed. For a bunched beam consisting of thousands of particles, the envelope equation is 
considered to understand the collective behavior of the beam. Then using PDT code to simulate the transverse dynamics 
of the tracking particles with the increase of energy, and to match FODO-like lattice periodic beam channel. 

Analyze the transverse force on the individual particles, and substitute the radial electric field and angular magnetic 


field formulas derived in formula (3). The dynamic equations can be obtained: 


dp, 9 (= F aa) 
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ae as E (5) 
dz MoVz 

MoVz 


where q = e = 1.6 x 1071°C, the unit of p, is kg «m/s and: 


ðE, dA, doz |. 
== = TZ cos (Qo — $2(2) + wt) + sin (po — $a(2) + wt) a 
= = —wsin (o — $,(Z) + wt) 


According to the incompatibility theorem, in the process of acceleration longitudinally, the RF electric field presents 
the impact of defocusing to the beam transversely. To counteract this defocusing effect, external permanent magnetic- 
quadrupole lenses are needed. These lenses can be strategically placed to focus the beam and compensate for the 
defocusing caused by the RF electric field. By properly adjusting the strengths of these lenses, the beam can be kept in a 
tightly focused state, leading to improved transverse dynamics and reduced emittance growth. In the high-energy tank, 
the beam emittance is low enough that the FODO-like periodic structure can be adopted. However, it is important to note 
that in a proton Linac lattice, unlike a regular FODO lattice, the length of the accelerator structure grows proportionally 
to the relativistic beta of the accelerated particle. 

The transverse motion of a single proton is generally measured by four-dimensional parameters, including two 
position states x, y and two momentum states P,, Py or velocity states vx, vy or by the inclination x’ = v/v: = PsP}, y' = 
vy/v, = Py/P:. In addition, the Twiss parameter is generally used to describe the bunched beam envelope. The general 
formula for a phase space ellipse (an ellipse in the plane with axis x-x' and containing phase space particles for example) 


is as follows: 


f + cen) +( x i i (7) 
y ex/Vx y ExVx 


And the same is true in y-direction. Where y, a, p are Twiss parameters or so-called Courant-Snyder parameters, 
neither of which are independent parameters, but together constitute Courant-Snyder conditions:yB — a? = 

The distribution modes of the beam particles in phase space are Gaussian and Kapchinskij-Vladimirski (KV) 
distribution for typical numerical simulating. The more significant that need to be studied is the central and denser part 
of a beam while the outer parts of a beam with a very low density are usually not very important. The relation between 
the total emittance of n-dimensional uniform hyperellipsoid and rms emittance is that éo1a=€rms(n+2). When n=2, the 
points in the four-dimensional ellipsoid are uniformly distributed on all its two-dimensional projections, which is also 
known as KV distribution [21]. When x and y planes are both considered transversely, the four variables (x, x’, y, y^) are 


considered to obey KV distribution in the four-dimensional space, then there is étoai=4érms: 


F + cen) +( x ) P P + (Qay/yy yt +( y’ | = (8) 
V Ex/ Vx V ExYx véy/Yy VevYy 


In the PDT code, under the given Twiss parameters and emittances within limits over two transverse directions, take 
a required number (the input particle quantity) of groups of uniformly distributed pseudo-random numbers in a certain 
range for three parameters D1, ®2, ®3 E [0,27], which meet the condition: 
(cos ®, - cos ®,)* + (sin ®, - cos 3)? + (cos ®, - sin ®,)? + (sin®, - sin ®3)? = 1 (9) 
Then we can assign the values to the transverse state quantities (x, x’, y, y’) of particles: 


(ye / Vee 
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The input emittance is the single plane RMS emittance En. rms, which need to be converted into unnormalized 


= cos ®,-cos 中 ， 


= cos ®, -sin 中 ， 


= sin ®,-sin P; 


emittance Ems = En.rms/PoYo and then to the total emittance of hyperellipsoid sroral = 4Erms. According to the 
transverse dynamic formulas (5-6), calculate and track the motion state x, x’, y, y’ of each particle with PDT code, of 
which we use difference method to calculate the partial derivative of electric field. 

After tracking a bunch of particles with typical input parameters for a test, it is observed that the particles pass 
through a continuous periodic RF structure with an accelerating gradient of 30 MV/m. The resulting phase space 
distribution diagram, shown in Figure 2, indicates that the bunches have diverged noticeably. This divergence is evident 
in the x-x' phase space as larger displacements and velocity divergences, while in the y-y’ phase space, the beam transitions 
from a convergent state with a negative parameter a to a divergent state with a positive a. 
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Fig.2 (color online) The distribution of the particles in unnormalized transverse phase-space at the beginning (blue) and at the end (red) of 


the acceleration 


Excessive transverse drift distance will result in increased particle losses. To deal with it, it is essential to introduce 
elements with focusing capability. A FODO-like lattice is composed of two acceleration tanks and two PMQs in such 
structure in alternating arrangement, where each tank consists of several identical RF structure cells. To find the matching 
point of the Twiss parameters to the gradient of the PMQs, a process known as matching, it is necessary to ensure that 
the beam's transverse properties align with the focusing strength provided by the PMQs. 

The linearization of the change in transverse position and velocity inclination of the particle between the two 


positions can be expressed by a 2 X 2 transfer matrix M. Then the transfer matrix of the accelerating cavity can be 


obtained in this way. Assume that P, = x’: By, then we have the normalized motion state which is: 
-1 


1 0 
“2 | = M. ff ewes 
Pal =o Mlo 1| [pa] =m leal aD) 
B2Y2 Biy 
For periodic FODO structure, its transfer matrix can be expressed by Twiss parameter: 
M= - + asin Bsinpg | \M*| =1 (12) 
—ysinu cosy — asinu 


In the case of calculating the transfer matrices for the accelerating cavity, the calculation is based on the state quantity 
of the particles. These transfer matrices are independent of the motion of individual particles and the presence of the 
PMQs. Instead, they solely depend on the structure of the accelerating cavity. Therefore, once the structure of each tank 
is determined, the corresponding transmission matrix, denoted as M*, can be calculated using equation (12). 

Regarding the PMQs, its characteristics dictate that it cannot simultaneously focus or defocus in both the x and y 


directions. The transfer matrix associated with the PMQs differ between the two transverse directions: 


sinVK1 sinhvK1 
VKL 一 一 hyk 一 

Rfocus = VK | Rdefocus = VK (13) 
一 VKsinVKI cosVKI VKsinhVKI coshVKI 


where K = qG/mcBy = q(B/r)/mcPy, G is the gradient of the quadrupole lens (in T/m), B is magnetic field intensity 
at pole tips, / is the length of the lens. 
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Fig.3 (color online) Solving procedure to match the beam in a FODO-like lattice 


Fig.3 shows the matching process. A FODO-like lattice transfer matrix is changed by adjusting the parameters G 
and G2 of two quadrupole lenses with fixed parameters of the acceleration cavities. The Twiss parameters, for depicting 
the envelope of the phase-space ellipse, obtained by matching the import beam to export in this lattice after normalization. 
The next steps are calculated after un-normalization. 

In the optimization process, numerically, there will be such a situation that the trace of the matrix |Mj, + M3,| > 2, 
in which the matching solution or a matched beam condition does not exist. Hence, there’s need to adjust the gradients 
G1 and G2 of the two quadrupoles lens so that the x direction meets |Mj, + M3.| < 2, and the same is true in y direction. 


The code seeks results automatically after the end of the longitudinal dynamics design, allowing determination of optimal 


Twiss parameters of the phase space ellipse and gradients of the magnetic-quadrupole lens which leads to best 
transmission. Ultimately, in a testing FODO-like lattice structure, the normalized envelope of phase-space ellipse in the 


two directions is shown in Fig.4. It reveals that the normalized ellipse almost matches perfectly. 
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Fig.4 (color online) The distribution of the particles in normalized transverse phase-space at the beginning (blue) and at the end (red) of an 


optimal matched FODO-like lattice solved by PDT 


3 Benchmark 


Here a coupled cavity linac (CCL) structure is considered for the benchmark against PDT and Tracewin, which 
demands as an input the standing wave field map. Nevertheless, standing wave is actually a special form of travelling 
wave. CST is used to generate the electromagnetic field map of CCL structure for a single cell. 

In this scenario, a group of 10,000 protons with an incident kinetic energy of 85.8 MeV undergo transportation 
through a FODO lattice. The lattice consists of two 7-cell CCL accelerating tanks responsible for boosting energy and 
two PMQs for focusing. The system operates at a frequency of 2.97863 GHz, and the dimensions of the acceleration 
structure can be found in Table 1. Given the same particle incident phase 117 deg and initial emittance of 0.3978 
TImm'mrad, the optimal matching values are obtained respectively by automatic optimization. The results demonstrate 
reliability and accuracy in terms of kinetic energy growth, transverse emittance, and the gradients of the PMQs. 


Table 1 The specific parameters of benchmark 


Emittance (Norm.) Output PMQs 
Acceleration structure length [m] 
[xmm'mrad] energy [MeV] gradients [T/m] 
PDT 0.398 89.506 445/445 CCL cell PMQs Drift Active/Total 
TraceWin 0.400 89.457 446/447 0.0202 0.02 0.01 0.2828/0.3622 


The phase space distribution of particles, both in terms of envelope and Twiss parameters, exhibits a good fit with 
the results obtained from both PDT and TraceWin, and is shown in Fig.5. Furthermore, the variation of the final energy 
of the particles with the input phase also exhibits a high level of agreement, with an error within 0.02%. These findings 


showcase the accuracy and effectiveness of the PDT code used in the study. 
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Fig.5 (color online) Normalized emittance and envelope of the phase-space ellipse (left) 


Kinetic energy grows with accelerating and final energy varies with input phase (right) 


4 Dynamics simulation of a 30MeV-100MeV multi-tanks BTW linac 


The CS-30 cyclotron, in operation by the Institute of Nuclear Science and Technology of Sichuan University, can 
accelerate protons, deuterons and a particles, the energy of extracted proton beams is about 30 MeV [22], is a good 
candidate for high gradient acceleration structure test. What we expect is to take proton beams to a higher energy to make 
it more serviceable in wider applications. The author has designed a S-band side-coupled cavity linac (SCL) for re- 
accelerating proton beam from 26 MeV up to 120 MeV, its shunt impedance ranged from 22.5 to 59.8 MQ/m, and the 
designed acceleration gradient is about 13 MV/m, however suffering hardship of the tuning of the cells [23]. The BTW 
exhibits superior performance in terms of shorter filling time, higher TT factor, sharper nose angle in the nose region, 
and approximate values of ZTT [24-25], contributes to make feasible a compact S-band BTW proton linac. 

Protons are decided to be injected into a BTW linac with an acceleration gradient of 30 MV/m, operating at f = 
2.856GHz. The linac is composed of 10 accelerating tanks, and each tank consists of multiple BTW cells of equal length. 
In the process of acceleration, there can be transverse divergence of particles due to the RF fields, and PDT code 


optimization algorithm is used to make transversely beam matching. 


4.1 BTW structural parameters tuning 


T 


Fig.6 (color online) Model diagram of BTW accelerating structure 


The geometry of the regular BTW cell structure is depicted in Fig.6. It is recognized that the thickness of the iris 
influences the shunt impedance and mechanical resistance of the structure remarkably. We can adjust the wave velocity 
by varying the cell length. In addition, group velocity is related to the radius of coupling holes. The group velocity is 
changed by changing the size of the coupling holes for travelling wave structure, which is located far away from the beam 


axis, and it is reasonable to assume that the field along axis is unaffected by the size of coupling holes. The RF frequency 
is related to the cell cavity diameter and the radial distance of coupling holes. Fig.7 illustrates the longitudinal field 
amplitude and phase of a two-cell constant gradient structure constructed from single cell calculation results. 
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Fig.7 (color online) The amplitude and phase of the electric field of the adjacent units under periodic cells 


The relation between the cell length Lee and the phase velocity py of the BTW cell, the phase advance o and 
frequency f of travelling wave is as follows: 

Oppc 

Leen = Daf 

nf 

f=2.856GHz and o=5z/6 of each cell are constant, the length Lee has a linear relation with phase velocity £p. The 


(14) 


initial energy of the particle is 30 MeV, where we have particles initial velocity 6o=0.247. When the acceleration of first 
tank is finished, the particles enter the next cavity for acceleration with velocity increasing significantly. Therefore, the 
phase velocity of the travelling wave in the next tank should be increased accordingly to make the it consistent with the 
particle velocity as far as possible and follow the rule between fp and £o in the meantime, so as to achieve the highest 
acceleration efficiency. 

The transit-time factor (TTF) is the ratio of the actual energy gain in an acceleration cavity to the energy gain 


obtained by the specific dc field, and represents the acceleration efficiency of the cavity [21]: 
L 
Í Az, real COS (az 一 Pz + go)dz 
L 
Í Az, real dz 


The TTF value is affected by the initial velocity, incident phase, the number of cells in each tank and travelling wave 


TTF = (15) 


phase velocity. The left figure in Fig.8 shows the change of the maximum TTF with the number of cells and phase velocity 
fp of travelling wave and the parameters of the acceleration cavity change correspondingly with it, enable finding the 
optimal value for number of the cells in a tank. On the other hand, it is hoped that each tank contains as many cells as 
possible to make the particle energy reach 100 MeV as soon as possible. A compromise is considered to select 17 BTW 
cells in each tank. 

When keeping the parameters of cavity constant, the ratio of the energy gain with respect to the maximum energy 
gain V/Vinax of the reference particle in a tank varies with the particle initial velocity fo and the number of cells, as shown 
in the Fig.8 on the right. In the graph, it is observed that when the phase velocity of wave is constant £,=0.257, the 
influence of the number of cells on the optimal incident velocity of particles is not significant. And the optimal incident 
velocity /o=0.247 of particles in the first tank can be obtained, corresponding to the incident energy of 30 MeV. 
Conversely, the optimal phase velocity Jp, of travelling wave corresponding to different incident velocities Jo can be 
obtained in this way. Through extensive analysis of experimental data, it is convincingly demonstrated that the optimal 
phase velocity £, should be 0.01 larger than the particle velocity fo to achieve maximum energy gain throughout the 


entire acceleration process. 
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Fig.8 (color online) The maximum TTF varies with the Number of cells and phase velocity £p (left) 
V/Vinax varies with the Number of cells and £o, holding fp constant (right) 


The electric field on the axis of cavity and the optimal cavity geometric dimension are calculated by CST. The 
frequency, denoted as f, was adjusted by modifying the outer radius of the cell cavity, Rou, while maintaining a constant 
ratio of tube length (Lise) to cell length (Leen) at 0.6. And ensured that f remained at 2.856 GHz. Subsequently, the PDT 
code was utilized to calculate the optimal incident phase and energy gain as the particles passed through a 17-cell tank. 
The obtained parameters, such as the velocity p, cell length Leen, and cavity nose length Lie, were then used as inputs 
for CST in the next iteration. This iterative process was automated and repeated until the particle energy reached 100 
MeV. 

At the conclusion of the structure tuning simulation, the particle energy was accelerated to approximately 102.7 
MeV. The specific parameters of the accelerating cavity can be found in Table 2. 


Table 2 The specific parameters in each tank 


Particle initial Phase Output energy Transit time Cell length 
™ velocity fo velocity py [MeV] factor (TTF) [mm] 
1 0.2470 0.2570 35.2910 0.884 11.8994 
2 0.2668 0.2768 41.0496 0.889 12.8290 
3 0.2865 0.2965 47.2542 0.894 13.7481 
4 0.3059 0.3159 53.8974 0.898 14.6554 
5 0.3251 0.3351 60.9820 0.902 15.5495 
6 0.3440 0.3540 68.5051 0.906 16.4290 
7 0.3626 0.3726 76.4541 0.909 17.2937 
8 0.3808 0.3908 84.8124 0.912 18.1417 
9 0.3987 0.4087 93.5835 0.914 18.9721 
10 0.4161 0.4261 102.7443 0.917 19.7841 


The optimal incident phase is obtained through the optimization algorithm in PDT code. The TTF of particles in 
each tank is above 0.88, which achieves the desired acceleration efficiency. The total length of BTW structure is 2.708 


m, which is designed for the main acceleration cavity. 
4.2 Beam Matching 


As far as the transverse divergency of the beam is concerned, the magnetic-quadrupole lenses are incorporated in. 
According to the design of the main acceleration section, the structure of the acceleration cavity is determined. PDT code 
is used to solve the magnetic field gradient of two quadrupole lenses with thickness of /=0.03m in the first lattice, and 
Twiss parameters and the gradient is obtained G;=G2=3847/m. In order to facilitate the experimental control and 


engineering manufacturing of PMQs, the other eight lens are set to be exactly the same as in the first lattice. The transverse 


ellipses can almost be matched after five complete FODO-like lattice configurations. Fig.10 shows the normalized phase- 
space ellipses at the beginning and the end of the linac and the envelope curves are successfully limited to an acceptable 
range. 
The normalized emittance, as shown in Fig.9, remained at about 0.1 2-mm- mrad throughout the accelerating process, 
with a very small emittance growth rate and a complete beam transmission. 
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Fig.9 (color online) Normalized transverse emittance variation 


Then the longitudinal multi-particle tracking simulation is carried out. The maximum longitudinal phase angle offset 
is assumed to be Ag = +5°and the initial maximum energy difference is AW = +0.1MeYV relative to the incident 
reference particle. And ideally, the particles are evenly distributed in the longitudinal phase space. After the overall 
accelerating, the phase space envelope diagram of the longitudinal bunched beam is shown in Fig.10. It is found that the 
ellipse of particles within a certain range of the phase and energy from reference particle, when the acceleration period 


has terminated the phase width occupied by the bunched beam compresses in the longitudinal direction, while the energy 
spread increases. 
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Fig.10 (color online) Configuration of the overall structure of linac (F and D are quadrupole lenses), with the envelope throughout the 
acceleration (Xmax, Ymax, A@ and AW) and the normalized phase-space ellipses at the beginning and the end of the acceleration (T is 


transverse while L is longitude representing on the same axes respectively) 


5 Conclusion 


This paper presents the development of a new particle dynamics tracking (PDT) code that optimizes the dimensions 
of the overall accelerating structure and establishes electromagnetic fields in different acceleration tanks according to the 
particle motion with CST and MATLAB co-simulation. It supports both standing wave and travelling wave 
configurations, making it a versatile tool for a variety of particle types. Additionally, the code allows for the inclusion of 
additional elements, such as PMQs or solenoid for transverse beam focusing. And a benchmark is successfully performed. 

The proposed high gradient proton linac, based on the BTW structure, is divided into 10 acceleration tanks with 
increasing length. Between the tanks, 10 PMQs are placed for transverse focusing of the bunched beam. The linac 
stretches a total length of approximately 3 m, and the overall configuration of the structure is depicted in Fig.10. The 
target average gradient for the BTW structure is 30 MV/m, and the overall average TTF of the acceleration cavity is 0.9, 


with the characteristics of full transmission and minimum emittance growth. 
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质子 行 波 直线 加 速 器 的 束 流动 力学 模拟 程序 设计 


HE, TER, HHL, ESE, PB 


K, FAR 


By 


(辐射 物理 及 技术 教育 部 重点 实验 室 ， 原 子 核 科学 技术 研究 


摘要 : ”和 荷 能 质子 束 在 癌症 治疗 、 工 业 发 展 和 科学 研究 等 方面 的 广泛 应 用 


所 ， 四 川 大 学 ; 成 都 610065) 
， 使 得 质子 加 速 器 的 研究 日 益 升 温 。 近 年 来 ， 


数 GHz 的 高 梯度 射频 驻 波 和 行 波 结构 的 发 展 为 紧凑 型 粒子 直线 加 速 器 的 研究 提供 


子 直线 加 速 器 设计 工具 相 比 ， 行 波 直线 加 速 器 的 仿真 设计 工具 非常 有 限 。 
跟踪 (PDT) 程序 来 模拟 粒子 的 状态 轨迹 ， 同 时 满足 对 整体 


加 速 链 各 个 元 件 的 自动 调谐 。 其 模块 化 结构 使 其 可 以 根据 仿 


了 新 的 思路 。 然 而 ， 与 丰富 的 驻 波 粒 
因此 开发 了 一 个 基于 行 波 结构 的 粒子 动力 学 


= 


动力 学 仿真 程序 进 


真 设计 需求 调整 添加 功能 ， 包 括 对 类 FODO 结构 束 流 匹配 以 及 寻求 能 量 调制 需要 的 最 优 起 点 等 。 与 TraceWin 直线 束 流 


行 了 比较 测试 ， 结 果 表 明 该 方法 在 能 量 增益 和 束 流 匹配 上 都 具有 良好 的 精度 。 


该 程序 已 应 用 于 一 个 


工作 频率 为 2.856 GHz 的 低 流 强 返 行 波 质子 直线 加 速 器 的 束 流动 力学 设 让 
段 和 10 个 四 极 永 磁 磁 铁 组 成 ， 加 速 梯度 为 30 MV/m， 具 有 
关键 词 : 质子 束 ; 行 波 ， 粒 子 动力 学 ， 束 流 匹 配 ; 直线 加 速 器 
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EF 中。 该 直线 加 速 器 全 长 约 3 米 ， 由 10 个 加 速 


完全 传输 和 低 发 射 度 增长 的 特点 。 


